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ABSTRACT 



This study compares three single variable search methods — 
Golden Section, cubic interpolation and quadratic interpolation. 
The SUMT nonlinear program was used for the comparison. The 
OPT subroutine which performs the single variable search in 
SUMT currently uses the Golden Section method. Two different 
OPT subroutines were written which implemented cubic inter- 
polation and quadratic interpolation. Seven test problems 
which contained 9-100 variables and 2-20 constraints were used. 
The comparison was made on computation time per single variable 
search for the three methods and the number of function evalua- 
tions per single variable search for the Golden Section and 
quadratic interpolation methods. 

A single variable search by Lasdon, Fox and Ratner and 
one by Fletcher and McCann were also discussed. 

The results showed that the quadratic interpolation was 
slightly faster than the other two methods and required fewer 
function evaluations per single variable search than the 
Golden Section method. Time per single variable search was 
approximately the same for the cubic interpolation and Golden 
Section methods. Cubic interpolation required fewer points to 
be evaluated than the other two methods but the need for 
gradient evaluations proved to be costly in terms of computation 
time per single variable search. 



4 



TABLE OF CONTENTS 



I. INTRODUCTION 8 

II. THE PENALTY FUNCTION METHOD — 10 

A. THE OPTIMIZATION PROBLEM 10 

B. THE SUMT COMPUTER PROGRAM 10 

III. SINGLE VARIABLE SEARCH METHODS 18 

A. GOLDEN SECTION SEARCH METHOD 19 

B. CUBIC INTERPOLATION METHOD 23 

C. QUADRATIC INTERPOLATION METHOD 28 

D. SECTIONS COMMON TO THE DIFFERENT METHODS 32 

E. OTHER SINGLE VARIABLE SEARCH METHODS 34 

1. Lasdon, Fox and Ratner Method 34 

2 . Fletcher-McCann Method 39 

IV. TEST PROBLEMS 44 

A. TEST PROBLEM ORIGIN 44 

B. TEST PROBLEM DESCRIPTION 44 

V. TEST RESULTS AND COMPARISONS 51 

A. PROGRAMMING AND TESTING PROCEDURE 51 

B. COMPARISON CRITERIA 52 

C. COMPARISONS 52 

VI. CONCLUSIONS AND SUGGESTIONS FOR FURTHER STUDY — 57 

A. CONCLUSIONS 57 

B. SUGGESTIONS FOR FURTHER STUDY 58 

APPENDIX: COMPUTER CODES AND GLOSSARIES 60 

LIST OF REFERENCES 75 

INITIAL DISTRIBUTION LIST 76 



5 



LIST OF TABLES 



I. Test Problem Results 53 

II. Normalized Test Problem Results 55 



6 



LIST OF FIGURES 



1. Behavior of the Penalty Terms 13 

2. SUMT Program Flow Diagram 16 

3. Golden Section OPT Subroutine 22 

4. Cubic Interpolation OPT Subroutine 27 

5. Quadratic Interpolation OPT Subroutine 31 



7 



I. 



INTRODUCTION 



The concept of optimization has grown to play a major 
role in the analysis of many complex decision and allocation 
problems. The quality of the analysis not only depends upon 
the skill and good judgement used in interpreting and 
modeling the problem but also upon the reliability and 
efficiency of the vehicle used for finding a solution to the 
problem. The solving of nonlinear programming problems has 
seen substantial development during the past twenty years 
and, no doubt, will increase in acceptance and popularity in 
the future. 

Within the Department of Defense several problems of 
significant importance, such as weapons targeting and alloca- 
tion, aircraft and power plant design and manpower planning, 
to name a few, have arisen which are nonlinear in nature. 

Some examples of industrial nonlinear programming problems 
include oil refining, curve fitting, inventory and logistics. 

Many methods exist for solving nonlinear programming 
problems. Some of these methods, called penalty function 
methods, are based on transforming a constrained minimization 
problem into a sequence of unconstrained minimization 
problems whose solutions approach the solution of the 
original constrained problem. 

In the penalty function method, as well as in most other 
nonlinear programming algorithms, a one dimensional 
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minimization search is used repeatedly in the solution 
process. In this paper several different one dimensional 
minimization search methods are compared in the context of 
penalty function algorithms . 

The different one dimensional minimization search methods 
(also called single variable searches) are compared by the 
amount of computation time and by the number of function, 
gradient and Hessian evaluations of the objective function 
and constraints needed to find a solution to the original 
constrained optimization problem. 

Section II of this paper explains the formulation of 
the penalty function method and where the single variable 
search is used in this method. The nonlinear computer 
program SUMT is used in comparing* the different single 
variable search methods and is also explained. 

Section III formulates the different single variable 
search methods — Golden Section, cubic interpolation and 
quadratic interpolation — and also discusses several other 
methods which were not programmed because of their 
complexity. 

Results and comparisons of the different methods are 
stated in Section IV with conclusions and suggestions for 
further study in Section V. The Appendix contains the 
computer codes for the three methods that were programmed. 
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II. THE PENALTY FUNCTION METHOD 



A. THE OPTIMIZATION PROBLEM 

The single variable search is a major part of finding 
an optimal solution to a constrained optimization problem 
by the penalty function method. Dr. W.C. My lander, an 
author of the SUMT computer code, estimated that forty percent 
of the computation time required to find an optimization 
problem solution using SUMT involved single variable searches. 

The constrained optimization problem to be solved is of 
the form: 

minimize: f(x) 

subject to: gj (x) ^0 j = l,...,m 

hj (x) = 0 j = m+l,...,m+p , 

where x is an n-dimensional vector, the functions f, g and h 
are continuous and may be either linear or nonlinear, and 
the set of values satisfying the inequality constraints has 
a non-empty interior, i.e. {x: g^ (x) > 0, j=l,...,m) . 

B. THE SUMT COMPUTER PROGRAM 

The computer program used to compare the different single 
variable search methods was SUMT - Version 4 (Sequential 
Unconstrained Minimization Technique for Nonlinear Programming) 
coded in FORTRAN IV by W.C. My lander, R.L. Holmes and G.P. 
McCormick for the Research Analysis Corporation [Ref. 1] . 
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It was written to experiment with different methods of 
solving nonlinear programming problems. 

A Guide to SUMT - Version 4 [Ref. 2] is the program manual 
which contains detailed information about the development of 
the code, the subroutines, options, input-output, limitations, 
user-supplied subroutines, data deck setup and a detailed 
example of the use of the program to solve a nonlinear 
programming problem. 

The method used by SUMT to solve the original constrained 
optimization problem is described in detail in Chapter 8 of 
the book on nonlinear programming by Fiacco and McCormick 
[Ref. 3] . Basically, SUMT solves the constrained problem 
by finding the solutions of a sequence of unconstrained 
problems which approach the solution of the original con- 
strained problem. 

Previous versions of SUMT used different penalty functions 
for approximating the solution to the constrained problem but 
the function currently used to sequentially solve the problem 
is of the form: 

m m+p 2 

P (x,r) = f ( x) - r I In g. (x) + S [h. (x)]Vr , 

j=l J j=m+l J 

where the parameter r determines the severity of the penalty 
and consequently how closely the unconstrained problem 
solution approaches the solution to the constrained problem. 

The penalty function of SUMT uses the mixed interior 
point-exterior point method described in detail in Chapter 4 
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m 

of Ref. 2. The term, r Z In g . (x) , involving the 

j=l 3 

inequality constraints is from the interior point method 

which is also known as the barrier method. As the x vector 

approaches the boundary of the infeasible region for these 

constraints, this term approaches infinity. The term, 
m+p 2 

E [h . (x) ] /r , involving the equality constraints is 
j=m+l ^ 

from the exterior point method. As the x vector moves 
farther away from the feasible region for these constraints, 
this term approaches infinity. The behavior of these two 
terms is shown graphically in Figure 1. How much of a 
penalty these terms add to the penalty function depends upon 
the value of the parameter r. 

The solution procedure requires minimization of P(x,r) 
over those x satisfying the conditions gj (x) > 0 , j=l,...,m 
for r = r ]_' r 2 '*** w here r^ < < . . . < r^. < ... <0 . Under 

mild conditions on the original NLP problem the minima of 
the sequential unconstrained problems approach the solution 
of the original constrained problem as r^ 0 . The condi- 
tions to be met for the method to solve the original problem 
are described in detail in the SUMT program manual [Ref. 2]. 

When the NLP problem is convex the SUMT program also 
generates a sequence of points that are feasible for an 
optimization problem called the dual of the original prob- 
lem. The. maximum function value of the dual feasible points 
and minimum function value of the unconstrained problem 
bracket the optimal solution of the constrained problem. 
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a. INTERIOR POINT METHOD INVOLVING INEQUALITY CONSTRAINTS 




b. EXTERIOR POINT METHOD INVOLVING EQUALITY CONSTRAINTS 
FIGURE 1. BEHAVIOR OF THE PENALTY TERMS 
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As the unconstrained problem is minimized and the number of 
dual feasible points found increases the difference between 
the two solutions decreases and the optimal solution is more 
accurately approximated. 

SUMT makes use of the theory derived in Fiacco and 
McCormick's book [Ref. 3, Chapter 5] that uses the Lagrange 
extrapolation technique to approximate the solution of the 
constrained problem when more than one minimum of the 
unconstrained problem has been found. By using these 
extrapolation approximations the solution converges faster 
and achieves a closer approximation to the actual solution 
than the minima of the unconstrained problems. 

The logic of the computer program taken from the SUMT 
manual [Ref. 2] follows. 

STEP 1 . 

Find an initial starting point x Q within the interior 
feasible region such that g^ (x) > 0 , j=l,...,m . If such 

an initial point is not given or is not feasible, the program 
uses the optimization method to find one. 

j 

STEP 2 . 

Determine r^, the initial value of r, which can either 
be an input parameter or a rule for choosing r^ can be 
specified. 

STEP 3 . 

Determine the minimum of the unconstrained penalty 
function for the current value of r. 
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STEP 4. 



Estimate a better solution using the extrapolation 
technique, if possible. 

STEP 5 . 

Computation is terminated if the stopping criteria are 
satisfied thus yielding an approximate solution to the 
original problem with accuracy depending upon the stringency 
of the stopping criteria. 

STEP 6 . 

If the stopping criteria are not satisfied select 

r k+l < r k * 

STEP 7 . 

Go to Step 3. 

A simplified flow diagram of the computer logic taken 
from the SUMT manual [Ref. 2, p. 5] is shown in Figure 2. 

The single variable search problem is involved in Step 3 
of the procedure. In order to minimize the unconstrained 
penalty function the program chooses a search direction in 
n-space. It is believed that by searching in this direction 
from the starting point of the subproblem the penalty function 
will initially decrease. A one dimensional search is then 
made along that direction until a local minimum is found. 

At the solution to this one dimensional search problem a 
new search direction is computed and a new one dimensional 
search is performed. The one dimensional searches are 
repeated, each time with a new search direction, until the 
unconstrained problem for a given r is solved. 
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FIGURE 2. SUMT Program Flow Diagram [Ref. 2, p. 5] 
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SUMT allows the user to choose one of three procedures 
to determine the search direction — by Newton's Method using 
first and second partial derivatives of the unconstrained 
penalty function, by using the negative of the gradient of 
the penalty function or by a variable metric method which is 
taken from an algorithm of the Fiacco and McCormick book 
[Ref. 3, pp. 170-175]. 

Finding a solution to the unconstrained n dimensional 
problem thus requires several single variable searches. 
Therefore, the more efficient the single variable search is, 
the faster the computation time will be. Given a direction 
of search, the SUMT program uses the subroutine OPT to 
perform the single variable search. In this paper, the 
presently used single variable search method in the OPT 
subroutine is compared to other single variable search methods. 
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III. SINGLE VARIABLE SEARCH METHODS 



Since a major portion of finding a solution to the 
constrained optimization problem involves single variable 
searches, the SUMT program was used to compare several 
single variable search methods. 

The current OPT subroutine of SUMT uses the Golden 
Section search method to perform the single variable search. 
Two different single variable search methods were programmed 
to replace the current OPT subroutine. The first method 
uses cubic interpolation once the local minimum is 
bracketed to accelerate convergence to that local minimum 
and the second method programmed uses quadratic interpolation 
to accelerate convergence once three feasible points that 
bracket the local minimum are found. Single variable search 
methods by Lasdon, Fox and Ratner [Ref. 4] and Fletcher and 
McCann [Ref. 5] are also discussed. 

The object of the single variable search is to find a 
scalar 9 called the step length such that x* = x q + 9s 
where x* is the x vector that corresponds to the local 
minimum of the penalty function along the search direction s 
from an initial starting point x Q . For ease of notation 
penalty function values are shown as a function of step 
length 9 in each single variable search since the initial 
starting point x q , the search direction s and the parameter 
r are constant, i.e. P(9) = P(x Q +9*s,r). 
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The optimal step length could possibly be found by 
moving along the search direction in a random manner. 
However, this would not be a very logical way to attack the 
problem and would lead to many needless evaluations of the 
penalty function which would increase computation time. 

The OPT subroutine is supplied with a search direction 
and an initial point from which to search. It must call 
other subroutines of the SUMT program for functional, 
gradient or Hessian evaluations and is expected to return 
with a penalty function value and corresponding x vector 
which is the local minimum for that single variable search. 
It is assumed that the penalty function has a finite local 
minimum along the given search direction. 

The following subsections formulate and explain the 
different single variable search methods. 

A. GOLDEN SECTION SEARCH METHOD 

The Golden Section search method uses the limiting 
properties of the Fibonacci series [Ref. 6] . Fiacco and 
McCormick's book [Ref. 3] define the steps of the procedure 
used in the Golden Section method. 

STEP 1 . 

First an upper bound step length 0^ is obtained with the 
lower bound step length 0 L initially equal to zero. 9y is 
obtained by evaluating the penalty function at successive 
step lengths which are in the limiting Fibonacci ratio. 
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(l + /5~)/2 ~ 1.618 , that is, 0.. = ^ (1.618)^ , where 

i=0 

k is the smallest integer j such that 

3 Q j ~ 1 Q 

P [ Z (1.618) ] > P [ Z (1.618) *] 

1=0 1=0 

j “ 2 l 

with 0 T updated as Z (1.618) . 

L £=0 

STEP 2 . 

The interval bounding 0*, the optimal step length, is 
reduced by computing two other step lengths within the 
interval (0^0^. Their corresponding values are: 

e a = e L + O.382 (0 U - 0 L ) 

0^ = 0^ + 0.618(0^-0^) 

STEP 3 . 

The penalty function values of the two interior step 
lengths, 0 & and 0^, are then compared. 

STEP 4 . 

If P(0 ) < P(0,), then the optimal step length which 

3 JD 

corresponds to the local minimum of the single variable 
search should be in the interval (6-^/9^) if the penalty 
function is unimodal. Because of the fact that 
0.382/0.618 » 0.618, reassign 0 y = 0^, 0^ = 0 a and calculate 
a new 0 as computed in Step 2. Return to Step 3. 

CL 
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STEP 5. 



If P(9 a ) > PfQjj)/ then the optimal step size should be 

in the interval (6 .0 TT ). Reassign 0 T = 9 , 0 = 0. and 

a. u u a a D 

calculate a new 0^ as computed in Step 2. Return to Step 3. 
STEP 6 . 

If the penalty function values at both interior step 

lengths are equal, reassign 0 T = 0 , 0 rT = 0, and calculate 

ii a. u o 

a new 0 a and 0^ as computed before by returning to Step 2. 
STEP 7 . 

When the stopping criteria is satisfied as explained 
in subsection III.D, 0* is approximated by (0 L + 0 y )/2 
yielding an x vector x* = x Q + 0*S which corresponds to an 
approximated local minimum P(0*) of the single variable 
search. 

A flow diagram of the Golden Section method presently 

used in the OPT subroutine is shown in Figure 3. Certain 

parts of the computer code that were the same in each of 

the OPT subroutines are discussed in subsection III.D. 

The Golden Section OPT subroutine initially updates 0 T 

as it finds more feasible points by increasing the step 

length until it finds a point which is infeasible or where 

the penalty function value is larger than the penalty 

function value at the previous feasible point. If it finds 

an infeasible step length, the penalty function is assigned 

3 6 

a very high value (10 ) , the next to the last feasible 

step length is assigned 0 L and the last feasible step length 
assigned as the lower interior step length 9 . If the first 

cl 
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FIGURE 3. Golden Section OPT Subroutine 
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step length from the initial point is infeasible, the lower 
interior step length 9 a = 0.382 0 , since only one feasible 
step length (9 = 0) has been found. OPT then uses 9„ and 

9 a to compute the upper interior step length 
0k = 9 a + 0.382 (9 d - 9 a ) . The optimal step length is 
bracketed and the subroutine continues to reduce the 
interval of uncertainty by comparing penalty function values 
of the interior step lengths until it finds two values or 
an interval that satisfies the stopping criteria. 

In finding the minimum of each single variable search, 
the Golden Section method only requires penalty function 
evaluations which are compared to show where the penalty 
function decreases to a local minimum and then increases to 
infinity as the infeasible region boundary is approached. 

The factor by which the interval of uncertainty is reduced 
remains constant in this method. Therefore, any information 
about the behavior of the penalty function other than the 
fact that the optimal step length is somewhere within the 
interval of uncertainty is disregarded. 

B. CUBIC INTERPOLATION METHOD 

The cubic interpolation method was programmed for 
comparison with the Golden Section and quadratic interpola- 
tion methods. Basically, this method approximates the 
penalty function by a cubic fit and then finds the minimum 
of the cubic equation to approximate the optimal step 
length of the single variable search. 
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In order to approximate the penalty function by a cubic 
fit in the single variable search, two feasible step lengths 
which bracket the local minimum along with their directional 
derivative values are needed. This method actually consists 
of two distinct parts — the bracketing section to find the 
two feasible points and the cubic interpolation section to 
find a minimum of the cubic function. The single variable 
search method is supplied with an initial starting point x Q , 
the penalty function value at x q [P(0)], the gradient of 
the penalty function at x Q [VP(0)] and a search direction s. 

The steps of the cubic interpolation method follow. 

STEP 1 . 

Since the directional derivative at the starting point 
P'(0) is negative, by finding a step length that has a 
corresponding positive directional derivative, the local 
minimum along the given direction of search will be bracketed 
if the penalty function is unimodal. 

The step length is increased to find an upper step 
length 9 y until an n is found such that 

n i 

P' [ ( Z 2 1 ) - 1] > 0 . 

i=l 

If, while initially increasing the step length, an infeasible 
point is encountered, the increment by which the last step 
length was increased is halved and subtracted from the 
current step length. This is continued until a step length 
that is feasible is found. If its directional derivative is 
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negative, the increment is increased by one half the 
increment length and added to the current step length. 

This is continued until a step length with a positive 
directional derivative is found. 

The lower step length 0 T is always the last feasible 

1j 

step length encountered with a negative directional 
derivative . 

STEP 2 . 

Once two step lengths have been found that bracket the 
local minimum along the search direction, the cubic equation 
is used to compute a step length which corresponds to the 
minimum of the cubic fit. This step length is found by the 
equation: 



P’<e 0 ) + u 2 - d 1 

9 - e w - (By - 0 L ) [pi _ pt + 2Uj 



where 



P(6 - P(6 ) 

U ! = P *< 0 L ) + P,( V " 3 [ 0 L - 0 D " ] 



U 0 = [U, 



- p 



(0 L )P 



(0u } ] 



1/2 



STEP 3 . 

p' (0) is evaluated and if it is negative then reassign 

0 = 0 and the left side is discarded. If P'(0) is positive 

Li 

then reassign 0y =0 and the right side is discarded. 
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STEP 4. 



If the stopping criteria is satisfied, the two penalty 
function values, P(0) and P(0 L ) or P(0) and P ( 0y) , depending 
upon P'(0), are compared and the smallest one is returned 
as the local minimum P ( 0 * ) of the single variable search 
with its corresponding x vector, x Q + 0*s. 

STEP 5 . 

If the stopping criteria is not satisfied return to 
Step 2. 

A flow diagram of the cubic interpolation OPT subroutine 
is shown in Figure 4. The parts of this subroutine that are 
the same in the other OPT subroutines are discussed in 
subsection III.D. 

The cubic interpolation method uses more information at 
each feasible step length to find the local minimum of the 
single variable search than the Golden Section or quadratic 
interpolation methods. It requires a gradient evaluation 
for the penalty function at each feasible step length in 
order to compute the directional derivative. A gradient 
evaluation may require many function evaluations depending 
upon the number of variables and constraints in the original 
constrained problem. However, with this extra information 
the local minimum should be found by evaluating fewer step 
lengths in the cubic interpolation section of this method. 

As the value of the parameter r decreases with each 
unconstrained problem, the penalty function rises more 
abruptly as the infeasible region boundary is approached. 
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FIGURE 4. Cubic Interpolation OPT Subroutine 
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This makes it more difficult to find a step length with a 
positive directional derivative which is required to bracket 
the minimum. 

C. QUADRATIC INTERPOLATION METHOD 

A quadratic interpolation single variable search was also 
programmed for comparison with the two previously discussed 
methods . 

This method requires three feasible step lengths and 

their penalty function values which obey the relationship 

P(0 f ) > P(0„) < P ( 0 TT ) where 0 T < 0 SJI < 6„ . With this 
L M U L M U 

information a quadratic fit is made with its minimum 
approximating the minimum of the penalty function for the 
given search direction. The penalty function is evaluated 
at the interpolated step length which minimizes the quadratic 
function and compared to P ( 0^) allowing the interval of 
uncertainty to be decreased. The process is repeated until 
the minimum of the quadratic fit approximates the local 
minimum of the penalty function or the interval of uncertainty 
becomes small enough to satisfy the stopping criteria. 

This method also consists of two sections — the bracketing 
section and the quadratic interpolation section. The single 
variable search is supplied with an initial x vector x q , its 
penalty function value P(0) and a search direction s. 

The steps of the quadratic interpolation search method 
follow. 
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STEP 1. 



The step length is increased from the initial starting 
point to find an upper step length ©y until an n is found 
such that 

n . n-1 

P[ ( Z 2; - 1] > P[ ( Z 2 1 ) - 1] . 

i=l i=l 

If, while initially increasing the step length, an infeasible 
step length is encountered, the increment by which the last 
step length was increased is halved and subtracted from the 
current step length. This is continued until a feasible step 
length is found. If its functional value' is less than P(0 L ), 
the increment is increased by one half the increment length. 
This is continued until a -feasible step length is found that 
has a larger penalty function value than the previous step 
length. Before this is repeated, reassign 0 L = © M and 
©M = Qy/ always ensuring that P(8 L ) > P(0 M ) • 

If only one feasible step length is found in obtaining 
the upper step length, where in this case P(0^) > P(9^) 

(which indicates the minimum is bracketed) , the last 
increment is halved and subtracted from the upper step 
length which yields the interior step length 0 m . 

STEP 2 . 

Now that the local minimum is bracketed, the quadratic 
function is used to find a step length where the derivative 
of the quadratic fit vanishes. 
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i, j — 1,2,3 where 1 =* L, 2 = M, and 3 = U. 

STEP 3 . 

If 0 > 0^, then the left side is discarded and reassign 
= 9 m and 0^=0. If 0 < 0 M , then the right side is 

discarded and reassign 0 rT = 0„ and 0 W = 0. 

U M M 

STEP 4 . 

P(0 M ) is evaluated and if the stopping criteria. are 
satisfied, the smallest of P(0 L ), P(© M ) and P(9 U ) is returned 
as P(0*) with the corresponding x vector, x q + 0*s. 

STEP 5 . 

If the stopping criteria are not met, return to Step 2. 

A- flow diagram of the quadratic interpolation method OPT 
subroutine is shown in Figure 5 . 

The quadratic interpolation method only requires function 
evaluations to perform the single variable search. It uses 
the knowledge of three feasible step lengths and their 
penalty function values to approximate where the local 
minimum of the single variable search is located. However, 
it must find an upper feasible step length that is within the 
interval where the penalty function is increasing to infinity 
and also this feasible step length must have a penalty 
function value that is greater than that of the interior 
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FIGURE 5. Quadratic Interpolation OPT Subroutine 
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feasible step length. This causes the interval containing 
possible feasible upper step lengths to be reduced even more. 
As the parameter r decreases for each unconstrained sub- 
problem, the penalty function increases more abruptly from 
the minimum as the infeasible region boundary is approached. 
This also increases the difficulty in bracketing the minimum. 

Once the three feasible step lengths are found, the 
quadratic interpolation can converge to the local minimum 
thereby solving the single variable search. 

D. SECTIONS COMMON TO THE DIFFERENT METHODS 

Some parts of the three different OPT subroutines were 
made similar so that they could be compared under the same 
conditions . 

Essentially the same stopping criteria was used in each 

subroutine. When the ratios of the two x vectors that 

bracketed the minimum or ratios of their penalty function 

-7 

values were within 10 of one, the subroutine returned to 
the program with a local minimum and a corresponding x vector 
as the solution to the single variable search. These criteria 
were checked in the Golden Section method whenever two 
interior step lengths had been computed and in the inter- 
polation methods both after the minimum had been bracketed 
and also after each interpolated step length had been 
computed . 

In the three different methods, the subroutines also 
contained counters which would cause the search no stop 
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after ten interpolations in the cubic interpolation method, 
twenty interpolations in the quadratic interpolation method 
or twenty-five feasible points in the Golden Section method. 
In the test problems that were run these counters were never 
reached, thereby, never satisfying this stopping criterion. 

In the cubic interpolation method another stopping 
criterion was used along with the ratio tests. When the 
cosine of the angle between the gradient of the penalty 

-7 

function and the search direction vector was less than 10 
the subroutine returned with a local minimum and x vector. 

It was stated earlier that it was hoped that the penalty 

function would initially decrease in the given search 

direction from the initial starting point. In the three 

different methods, if the ratio of any element of the 

current x vector and the corresponding element in the initial 

-7 

x vector came within 10 of one, the single variable search 
was restarted in the opposite search direction. In the 
interpolation methods if after one hundred tries the minimum 
was not bracketed or the step length had been decreased more 
than twenty consecutive times, the search direction was also 
reversed and the single variable search was restarted. 

The three different OPT subroutines called the subroutine 
EVALU for penalty function evaluations . EVALU would return 
back a variable after each call up which would show if the 
step length was feasible. The cubic interpolation method 
called the subroutine GRAD for the gradient of the penalty 
function when it was needed to determine the directional 
derivative. 
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E. OTHER SINGLE VARIABLE SEARCH METHODS 

Two other single variable search methods were also 
looked at but were not programmed because of their complexity 
and their need for changing other portions of the SUMT 
program than just the OPT subroutine. 

1. Las don , Fox, and Ratner Method 

The single variable search developed by Lasdon, Fox 
and Ratner [Ref. 4] uses a penalty function that doesn't 
contain equality constraints and the penalty term for the 
inequalities is also different than the SUMT penalty function. 
The unconstrained minimization problem is: 

... m 1 
minimize P(x,r) = F(x) + r Z - — t — r* 

i=l G i U) 

where F is the objective function and the G^'s are the 
inequality constraints. This model could be modified to 
handle equality constraints and a penalty function like that 
used in SUMT. 

The single variable search method is in three distinct 
stages — linear approximation, quadratic approximation and 
cubic interpolation. Linear approximations are made for the 
objective function and each constraint using F(0), F'(0), 

Gj ( 0) and Gj'(0). A step length (0^) which minimizes the 
penalty function consisting of the linear approximations is 
then calculated. The information used to make the linear 
approximations is then used along with F(0^) , F' (0^) , Gj(0^) 
and G j ' ( 0 -l ) to make quadratic approximations of the 
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objective function and constraints. The penalty function 
consisting of quadratic approximations is then minimized 
by computing 02* If the directional derivative of the 
original penalty function is positive at 0 2 > cubic inter- 
polation is applied yielding 0*, the solution to the single 
variable search. If the stopping criterion is satisfied 
during any stage, the single variable search is terminated 
yielding a local minimum and a solution to the search. 

The following steps outline the Lasdon, Fox and 
Ratner procedure. 

STAGE 1 . 

The objective function and inequality constraints 
of the original constrained problem are approximated by 
using the given values at the initial starting point. 

f 1 (9) = F(0) + F' (0) 

gij (0) = Gj ( 0) + Gj'(0) j = l,...,m. 

The approximating function for P(0) is 

m x 

Pl (6) = fl (e) + r ^ 

The smallest positive zero of the linear approximations of 
the inequality constraints is found by taking the smallest 
- Gj ( 0 ) /G j ' ( 0 ) over all j. This value is designated the 
upper step length 0 y . 

The step length corresponding to the minimum of the 
approximated penalty function is found by doing four or five 
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Newton's method iterations using i 9 y 



as the starting point. 



0 



1 




P l' ^2 

Pi" (i e D ) 



The result of this stage is a step length that 
approximates the optimal step length of the single variable 
search. If the stopping criterion, explained following 
Stage 3, is not satisfied proceed to Stage 2. 

STAGE 2 . 

The same procedure is performed in this stage as in 
Stage 1 with the exception that quadratic approximations to 
the objective function and constraints are made instead of 
linear approximations. The second stage is entered with 
values ’of F(0^) and Gj (0^) along with the information used 
to make the linear approximations. The quadratic approxima- 
tions to the objective function and constraints are 



where 



f 2 ( 0 ) = a9 2 + b© + c 



F(0 1 ) - b0 1 - c 




b = F* (0) , c = F ( 0) 



and 



where 



g (0) = d.0 2 + e. + f . j = 1/ ... /in 
^ 3 J J J 

G i <e 1 ) - e j6l - fj 

d . = — 1 « 

3 

ej = Gj ' ( 0) , f. = G.(0) . 
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If Gj(0^) is infeasible for some j, then the smallest 
positive root over all j of 



-e . ± /e . 2 - 4d . f . 

— 3 1 I_1 

2dj 



9 



2U 



is used as the starting step length for using Newton's 
method to minimize the approximated penalty function 



m 

p 9 (0) = f, (0) + r E 
2 i=l 




The step length which minimizes P 2 O) should be 
found after four or five iterations of Newton's method. 



0 



2 



P2’ (9 1> 

W 



If 0^ is infeasible, use j 0 2u in Newton's method. 
If the stopping criterion is not satisfied after finding 
©2 continue to Stage 3. 

STAGE 3 . 

Direct cubic interpolation of the penalty function 
is used in this stage. If P'(02) is positive, the minimum 
is bracketed and cubic interpolation can be performed to 
find a step length which corresponds to the minimum of the 
penalty function. If the minimum has not been bracketed, 
that is, if P' ( 0 2 ) < 0, a bracketing procedure like the one 

used in the cubic interpolation method of subsection II I. B 
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must be performed to find a step length that has a positive 
directional derivative. Once this step length is found, 
cubic interpolation is used to find a step length which 
minimizes the cubic function. If the new step length does 
not satisfy the stopping criteria, additional cubic inter- 
polations are made using the two points which most closely 
bracket the minimum of the penalty function. 

When two penalty function values are nearly equal 
or when the interval between step lengths that bracket the 
minimum becomes very small, an optimal step length is deter- 
mined by 

* = P'(e a )e b - p-(9 b )e a 

P' (0 ) - P' (0. ) 
a d 

where 0 is the lower step length and 0, is the upper step 

3 l O 

length (0 < 0* < 0^) . 

The stopping criterion used in this method is the 

I s^VP I 

— 1 1 , where <j> is the angle 

I I S | | | I VP I I 

between the search direction vector and gradient of the 
penalty function. When the cosine of the angle is less than 

— O 

10 , the stopping criterion is satisfied. 

The authors found that in using this method on test 
problems, the stopping criterion was often satisfied at the 
end of Stage 2. In their comparisons of this method with a 
cubic interpolation method, they found that it performed 
the single variable search much more efficiently [Ref. 4, 



cosine test | cos <J> | 



This method requires gradients of the objective 
function and constraints. In the SUMT program these values 
are not stored so implementing this method would require 
more storage or more gradient evaluations along with 
changes in other subroutines and/or addition of new 
subroutines . 

2 . Fletcher-McCann Method 

The Fletcher-McCann method [Ref. 5, pp. 210-212] 
uses an approximation of the original unconstrained penalty 
function to find the local minimum of the single variable 
search. The penalty function approximation used is of the 
form 



T ( 0 ) 



a + b0 + c0 



2 




/ 



where the parameters a, b, c, d and 0^ are determined by 
using objective function and constraint information at two 
feasible step lengths. The authors felt that an approximation 
of this type which goes to infinity at the barrier 0 = 0^ 
would fit the penalty function better than a polynomial 
approximation which goes to infinity only as 9+®. 0 y i s 
an estimate of the intersection of the search direction s 
with the boundary of the infeasible region. The approximated 
penalty function's minimum is found by applying Newton's 
method to find a step length which corresponds to the local 
minimum along the direction of search. 

Since the explanation of this method by the authors 
was rather sketchy, certain parts could be interpreted 
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differently and. only resolved by testing the method on 
different problems. An explanation of the notation for 
this method is first given. 



1 m+ P p 9 

F ( 0 ) = f ( 0 ) + ± E h . ( 0 ) s a + b0 + c0 2 

j=m+l 3 



m 

I 

j=l 



G(6) = -r I in 9 j (e) = 



F' (0) = 



p ? m+p m 

Vf (0) + ± E h. (0) s i Vh . (0) 
r j=m+l 3 3 

IUII 



~ b + c0 



G' (0) = - 



m s Vg . ( 0) 
E i 



j-i g : (e) 



(0 O - 6) 



The single variable search is supplied with an initial point, 
functional and gradient information at that point and a 
search direction. The steps of the method follow. 

STEP 1 . 

A feasible step length in the direction of search is 
needed along with the initial point in order to estimate the 
parameters a, b, c, d and 0^. A unit step length is initially 
taken in trying to find a feasible step length. If it is not 
feasible, the step length is progressively halved until a 
feasible step length is found. The lower and upper feasible 
step lengths are written as 0 Q and 0^, respectively. At the 
start of the search 0 q = 0. Function and gradient evaluations 
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are made at the second feasible step length in order to 
compute F^, G^, ' and G^' (shortened notation for F(0^), 

G ( 0 ^ ) , etc . ) . 

By using the information at the two feasible step 
lengths, simultaneous equations are solved to determine 
the parameters . The parameter values are 




b - V - 2ce o 



u 



t < e i-y 

v - G i 



The parameter 0^ is taken as the smallest value of two 
computed values which is greater than 0^. 

STEP 2 . 

Once the parameters are determined, the equation for 
the minimum of T which is 

T' = 0 = b + 2c0 + j 

(0 y - 0) 

is solved by Newton's method using the midpoint between the 
two feasible step lengths as the starting point. 
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0 



0 , - 0 . 






] 



where T" = 2c + \ . 

J (0U-6) 3 

STEP 3 . 

The values of VP(0) are computed and if the stopping 
criterion is satisfied, the single variable search is 
terminated with 0 as the optimal step length. 

STEP 4 . 

If the stopping criterion is not satisfied, a new 
step length is computed using Newton's method in Step 2. 

The authors of this method stated that different 
default actions would be necessary in cases where the inter- 
polation failed or was unsatisfactory. For example, if for 
the feasible step lengths, 0 q and 0^, a value of 0 y could 
not be determined, it would be necessary to use another 
feasible step length to determine where the intersection of 
the search direction and the boundary of the infeasible 
region was. In later stages of the minimization the authors 
said that this method sometimes failed and a quadratic 
interpolation would have to be used. 

Implementing this method in the SUMT program would 
also require more storage or more gradient evaluations 
along with modifications to other subroutines and/or additional 
subroutines . 
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No comparisons are made of the Lasdon, Fox and 
Ratner or Fletcher-McCann methods to the other three 
methods since they are only discussed and not programmed. 
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IV. TEST PROBLEMS 



A. TEST PROBLEM ORIGIN 

The problems used to compare the three different single 
variable search methods programmed as OPT subroutines were 
taken from problems used in a thesis by Lt. J. Waterman, USN, 
which compared three different nonlinear programming codes 
[Ref. 7], The SUMT program and problem data decks were the 
same as those used by Waterman except for the Golden Section 
OPT subroutine being replaced by the cubic and quadratic 
interpolation methods . 

The structure and degree of difficulty among the seven 
test problems are quite varied and represent a sample of 
some real world problems. The number of variables and 
constraints ranged from 9 to 100 and 2 to 20 respectively. 

The problems contain combinations of linear, nonlinear, 
equality and inequality constraints. 

The following subsection contains the descriptions of 
the test problems as presented in Waterman's thesis. 

Constants in each problem are represented by a, b, c, d, e, L, 
m, u and s unless otherwise specified. 

B. TEST PROBLEM DESCRIPTION 

Problem 1 was an example of determining the chemical 
composition of a complex mixture under conditions of 
chemical equilibrium. It contained 45 independent variables 
and 16 linear equality constraints. 
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minimize f(x) 



+ In 



x 



^ t ^ x ik^ C nk 
k=l j=l 3 * 3K 



jk 



Z x ik 
j=l 3 * 



-) ] 



subject to 



h ± (x) 



7 

Z ( Z E., k x ,) 
k=l j=l 1:|JC 1 3 k 



b i = 0, 



1 1 , . . . , 16 



x j 3 !/•••/ h k k 1 f • • • / 7 . 



Problem 2 was formulated by the Shell Development 
Company. It consisted of 15 variables and 5 nonlinear equality 
constraints . 



10 5 5 5 3 

maximize f(x)= Z b.x.- Z Z yz-2 Z d.z 

i=l 1 1 j=l i=l j=l 3 



where y = o ij x (10+i) and z = x (10+j) 



10 



subject 



to g.(x) = 2 Z y + 3 d . z + e. - Z a. .x. 0 

3 i=l 3 3 i=l 13 1 



x^ > 0 i - 1, . . o , 15 



Problem 3 was to maximize the area of a hexagon in which 
the maximum diameter was unity. The problem had 9 independent 
variables, 13 nonlinear inequality constraints and a lower 
bound of 0 for Xg . 
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Maximize: f(x) - . 5 (x^ - x 2 x 3 + x 3 x g - XgXg + XgXg - XgX ? ) 

subject to: 



2 2 

1 - x 3 - x 4 ■> 0 



1 - x 4 1 0 

i 2 2 

1 - x c - x- >0 
b 6 — 

1 - (x x - x 5 ) 2 - (x 2 -x g ) 2 > 0 
2 2 

1 - ^ - x 5 ) - (x 2 - Xg) 2: 0 

1 - (x x - x ? ) 2 - (x 2 - Xg) 2 _> 0 

1 - (x 3 - Xg) 2 - (x 4 - Xg) 2 _> 0 

1 - (x 3 - x 7 ) 2 - (x 4 - Xg) 2 2: 0 

1 ~ X ? 2 - (Xg - Xg) 2 2l 0 



x l x 4 ~ X 2 X 3 — ® 



X 3 X 9 1 0 



~ X 5 X 9 - 0 



X-.X-, - x-x,, > 0 
5 8 6 7 — 



x 9 1 0 • 



Problem 4 was probably the most difficult of the test 
problems. It included a linear objective function, 24 
variables, 12 nonlinear equality, 2 linear equality and 
6 non-linear inequality constraints. The variables had 
to also be zero or positive. 
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minimize: 



24 

f (x) = Z a . x. 

i=l 1 1 



subject to: 



h j .(x) = 


X 


( i+12) 


C i X i 




24 x. 

b Z — 1 

< i+12 > j=13 b j 


O 

tr 

H- 


X . 
_1 
b. 

3 


h 13 (x) 


24 
= Z 
i=l 


O 

II 

i— 1 
1 

X 






h 14 (x) 


12 
= Z 
i=l 


x. 24 

+ f Z 

d i i-13 


et ~ 1 - 671 
1 


= 0 




where 


f = 142.224 






U 


(x) = • 


" (x i + X (i+12) 


) 

+ & >0 


, i 


n (i+14) 


24 


I ^ 


9 



= 0 , i = 1, . . . , 12 



= 1,2,3 



j=l 



x . 
3 



h (i+14) (x) = 



(x (i+3) + X (i+15) } 
24 

Z x. 
j=l 3 



+ e^ >_ 0, i — 4,5,6 



x i il Or i = 1,. ..,24 
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Problem 5 was a weapon assignment problem with 100 
variables, a nonlinear objective function, 12 linear 
constraints and zero lower bounds for the variables. 



minimize : 



f (x) 



20 

E 

j=l 



Uj 



5 

n 

i=l 




1 ) 



subject to: 



5 

E x. . - b. > 0 j = 1,6,10,14,15,16,20 
i=l ~ 



20 

- E x.. + c. > 0 i = 1,...,5 

j=l 13 




i !,••«, 5, j — 1 



20 . 



Problem 6 was adapted from an inventory model where 
the x^, i = 1,...,50, represent the reorder quantity for 
50 inventory items and x^, i = 51,..., 100, represent the 
reorder points for the same 50 items. It contained 100 
variables, 1 linear and 1 nonlinear inequality constraint, 
and 50 lower bounds on the variables. 
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minimize : 



50 B. (x. + 50) 

f(x) = Z - 1 1 

i=l 



X. 

1 



where 



V X i + 50) =i (s i 2 +d i 2 ) *( i i) - 



s .d. d. 
i i A , lx 

2 ♦'i - ' ' 

1 



d i x (i+50) " m i * ^ (x) 



i_ e -x /2 



'2ir 



and $(r ) - f 4>(x) dx . 
t 



subject to: 



50 x. 

200,000 -^ c ihT + x (i + 50) 



- m i ) >_ 0 



50 L. 

300 E — > 0 
. , x. — 
i=l l 



X i > 0 , i = 1, . . . , 50 . 



Problem 7 was an entropy model. There were 46 population 
centers connected by a transportation network. Using a 
congestion cost function the model yields an equilibrium 
solution that identifies nodal populations as entropic 
functions of the total cost of the journey to work. The 
problem contained 46 variables, all of which have a zero 
lower bound, 1 nonlinear inequality and 1 linear equality 
constraint. 
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f (x) 



46 

E 

i=l 



minimize: 



X i x i 

500 (ln 500" * 



46 

subject to: 500 - E x- = 0 

i=l 1 



10000 - E c.y. + ad.v. > 0 
i=l 11 11 " 

where y. = x. + E x. , A(i) consists of all the 

jeA(i) 3 

arcs that converge directly and indirectly upon node i, 
and 

x^j>0 i = 1, . . . , 46 . 

The preceding problem descriptions are only meant to 
give a feeling of the kind of test problems used and their 
structure. Detailed information of how each problem was 
set up, the constant values, initial starting points and 
where the problems specifically originated is contained 
in Ref. 5. 
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V. TEST RESULTS AND COMPARISONS 



A. PROGRAMMING AND TESTING PROCEDURE 

The cubic interpolation method was first programmed 
as a program that did a single variable search when given 
a penalty function, its gradient, an initial starting point 
and a search direction. After it had been debugged it was 
converted into an OPT subroutine. After the subroutine 
was debugged using a couple of the less difficult test 
problems, OPT was converted into the quadratic interpolation 
method by modifying the bracketing procedure and changing 
the interpolation from a cubic function to a quadratic 
function . 

When the quadratic interpolation method had been debugged 
each problem was run through the SUMT program three times , 
using the three different search methods. As discussed 
earlier, the stopping criteria was made essentially the same 
in each of the OPT subroutines so that the only thing that 
varied in solving each test problem was the single variable 
search method. 

In the process of running the other test problems, minor 
debugging changes had to be made in the interpolation methods 
After all the changes had been incorporated, a final run 
of each problem with each method was made which produced the 
results used in the comparisons. Each method found the same 
solution to each problem to within six significant figures. 
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B. COMPARISON CRITERIA 



t 



Perhaps the best way to compare the different search 
methods is to compare the computation time required to find 
the solution to the problem. However, because of computer 
interactions, computation time may vary as much as twenty 

> 

five per cent when the same problem is run at two different 
times . Computation times were computed for each problem to 
see if a trend could be seen between the methods. 

Counters were inserted in the SUMT program so that the 
number of single variable searches performed and the number 
of functional and gradient evaluations needed to find the 
solution could be counted. 

Results of the test problem runs are shown in Table I. 

In running problem 7 with the interpolation methods it was 
necessary to change the factor by which the increment was 
multiplied for step length increase or reduction to 1.618 
vice 2 and 0.618 vice 0.5 respectively. This change was 
needed because a feasible starting point was found in the 
bracketing procedures which the SUMT program could not solve. 
By changing the factors to the same as those in the Golden 
Section method, the same feasible starting point was used 
to solve the problem. 

C. COMPARISONS 

Theoretically, each single variable search should reach 
the same solution for all three methods. However, because 
of the tolerances allowed in the stopping critera, the 
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TABLE I 



TEST PROBLEM RESULTS 



Problem 

1 

2 

3 

4 

5 

6 
7 

TIME = 
SVS = # 
F = # 
G = # 





Golden Section 


Cubic 


Quadrati 


TIME 


105.8 


105.7 


94.8 


SVS 


86 


89 


75 


F 


15487 


6171 


10149 


G 


1462 


7497 


1275 


TIME 


11.4 


9.6 


8.3 


SVS 


63 


53 


53 


F 


4354 


1574. 


2606 


G 


378 


1704 


318 


TIME 


4.8 


5.6 


4.0 


SVS 


39 


38 


40 


F 


8606 


3302 


5862 


G 


613 


3358 


628 


TIME 


395.3 


404.1 


396.1 


SVS 


163 


151 


174 


F 


44255 


13922 


30048 


G 


92.27 


97319 


98694 


TIME 


103.7 


97.9 


97.0 


SVS 


20 


19 


19 


F 


2683 


1250 


1991 


G 


296 


1415 


271 


TIME 


165.4 


145.9 


155.5 


SVS 


36 


34 


35 


F 


2320 


836 


978 


G 


104 


817 


107 


TIME 


78.8 


80.3 


67.8 


SVS 


20 


20 


17 


F 


865 


307 


397 


G 


1993 


2257 


1708 



Computation time 
single variable searches 
function evaluations 
gradient evaluations 
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solutions were slightly different at the conclusion of each 
single variable search. This meant that at the start of 
the next single variable search the search direction would 
be slightly different also. As several single variable 
searches were performed in each problem, the differences 
accumulated and this explains why for some of the problems 
it takes a different number of searches for each method to 
reach the same solution. This also explains the difference 
in the number of function and gradient evaluations required 
for each problem for the Golden Section and quadratic 
interpolation methods. 

Table II shows the normalized results for the time per 
single variable search, function evaluations per single 
variable search and gradient evaluations per single variable. 

As was expected, the number of gradient evaluations per 
single variable search is essentially the same for Golden 
Section and quadratic interpolation -methods . The only thing 
that can be compared between the three methods is the 
computation time per single variable search. Function 
evaluations per single variable search can only be used as 
a measure of effectiveness for Golden Section and quadratic 
interpolation since the cubic interpolation also requires 
gradient evaluations and fewer function evaluations. 

In looking at the times per single variable search the 
quadratic interpolation was faster than the Golden Section 
and cubic interpolation methods 5 out of 7 and 4 out of 7 
test problems respectively. Cubic interpolation was faster 
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TABLE II. NORMALIZED TEST PROBLEM RESULTS 



Problem 




Golden Section 


Cubic 


Quadratic 




TIME 


1.23 


1.19 


1.26 


1 


F 


180.1 


69.3 


135.3 




G 


17 


84.2 


17 




TIME 


.181 


.181 


.157 


2 


F 


69.1 


29.7 


49.2 




G 


6 


32.1 


6 




TIME 


.123 


.147 


.1 


3 


F 


220.7 


86.9 


146.6 




G 


15.7 


88.4 


15.7 




TIME 


2.42 


2.67 


2.28 


4 


F 


271.5 


92.2 


172.7 




G 


565.2 


644.5 


567.2 




TIME 


5.18 


5.15 


5.10 


5 


F 


64.4 


24.6 


27.9 




G 


2.8 


24.0 


3.0 




TIME 


4.59 


4.29 


4.44 


6 


F 


64.4 


24.6 


27.9 




G 


2.8 


24.0 


3.0 




TIME 


3.94 


4.02 


3.98 


7 


F 


43.3 


15.4 


23.4 




G 


99.7 


112.9 


100.4 



TIME = time/s ingle variable search 

= # function evaluations/single variable search 
= # gradient evaluations/single variable search 
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than Golden Section in 3 of the 7 problems with the same 
time for problem 2 . 

In comparing the number of function evaluations per 
single variable search, the quadratic interpolation method 
required fewer than the Golden Section on all seven problems. 

In most of the test problems, it required between thirty to 
forty per cent fewer function evaluations. The cubic inter- 
polation method required fewer function evaluations than the 
other two methods because it used more information about the 
penalty function at each feasible step length to find the 
optimal step length. However, it required more gradient 
evaluations which increased the computation time making it 
about the same as the Golden Section and slightly slower than 
the quadratic interpolation. 

With a larger sample of test problems more statistically 
sound comparisons could be made of the three different methods. 
Also samples of test problems that are similar in structure 
could be tested to show if one method worked better than the 
others on those type of problems. 

It should be emphasized that the results in this thesis 
are obtained from single variable searches on a very special 
type of function — the unconstrained penalty function in the 
mixed interior— exterior penalty function algorithm. No 
attempt should be made to generalize these results to other 
nonlinear programming methods which also use single variable 
searches, but on a significantly different class of functions. 
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VI. CONCLUSIONS AND SUGGESTIONS FOR FURTHER STUDY 



A. CONCLUSIONS 

In looking at the computation time required for each 
single variable search, the quadratic interpolation method 
was faster than the other two methods on the test problems 
that were run. The quadratic interpolation method uses 
information about the penalty function at each feasible step 
length to reduce the interval of uncertainty whereas the 
Golden Section reduces the interval of uncertainty by a 
constant factor. Because of this difference, quadratic 
interpolation proved to be slightly more efficient than 
Golden Section. The only drawback in the quadratic inter- 
polation method is the fact that bracketing the minimum can 
be quite hard to do. The cubic interpolation method reduced 
the interval of uncertainty in a more efficient manner than 
quadratic interpolation or Golden Section. However, having 
to make gradient evaluations to determine the directional 
derivative values proved to be very costly. The time per 
single variable search was about the same as the Golden 
Section search method. 

If the user of the single variable search method has 
prior knowledge about the complexity of the function or 
gradient evaluations he may prefer cubic interpolation over 
Golden Section or quadratic interpolation, or vice versa, 
since the cubic interpolation requires more gradient evalua- 
tions and fewer function evaluations. 
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B. SUGGESTIONS FOR FURTHER STUDY 

A combination of different methods made into a single 
variable search could be done with the three different 
methods that were programmed. This combined methods approach 
could use Golden Section until the minimum was bracketed, 
then use quadratic interpolation until the minimum 
was bracketed much closer and then finally use the cubic 
interpolation to home in on the local minimum of the search. 
The Lasdon, Fox and Ratner method uses three different stages 
that each find a feasible step length, with the final stage 
finding an optimal step length for the search. Another 
approximation of the penalty function like the Fletcher- 
McCann method could also be devised. 

Implementation of the Lasdon, Fox and Ratner or Fletcher- 
McCann method into the SUMT program with runs of the test 
problems would enable a comparison to be made not only to 
the three methods tested in this paper but also between the 
Lasdon, Fox and Ratner method and Fletcher-McCann method. 

Another possibility for study would be to fine tune the 
three different methods by adjusting the tolerances and 
termination conditions to see that if by allowing a less 
accurate determination of the minimum of the single variable 
search, the solution to the unconstrained problem could be 
found any faster. 

If a larger sample of test problems including more and 
varied types of nonlinear programming problems could be 
tested using the three different methods , a better comparison 
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could be made. The results may either show one method to 
be superior over the other methods in all of the problems 
or show each method suited best to a specific type of 
problem enabling the user to choose the method which best 
applies to the situation at hand. 
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APPENDIX 



DELX 

DELXO 

DOTT 

I 

ISW 

J 

KSW 

M 

MN 

N 

NSATIS 

NTCTR 

N401 

N404 

N405 

PREV3 



GLOSSARY FOR GOLDEN SECTION OPT SUBROUTINE [Ref. 2] 



The vector indicating the direction of move 
in one dimensional optimization. S 

The gradient vector of the penalty function 
VP. 



The inner product of the move vector and the 
negative gradient vector of the penalty 
function . -sT Vp 



Used as an index in DO loops. 

A switch showing whether the motion on a given 
vector failed to decrease the penalty function 
and the negative was tried. 

Used as an index in DO loops. 

Switch showing whether less than 25 feasible 
points were found on the direction vector. 

The number of inequality constraints 

Number of moves in search for a solution of 
a subproblem 

The number of variables 

Indicates whether constraints are satisfied. 

The number of the point on which the program 
is working. 

The number of points generated in attempting 
to find a point along the search direction. 

The number of infeasible points found on the 
direction vector. 

Switch showing that a feasible point on the 
direction vector could not be found and the 
negative was tried. 

Penalty function value at lower step length 
P(X3) . 
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PX1 



Penalty function value for left interior steplength 
P(X2). 

PO Current value of the penalty function P(X) 

Pi Penalty function value for upper steplength 

P(X1) 

P31 Penalty function value at initial starting point 

RJ Vector of current values of the constraints 

RJl Vector of previous value of the constraints 

X Vector of current value of the variables Xg+0s 

XX Temporary storage used for switching vector values 

XI X vector of upper step length X^+e^s 

X2 X vector of left interior step length W 

X3 X vector of lower step length W 
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SLERGUTINE OPT 
INFLICIT REAL*8 ( A-h ,0-Z ) 

C MARCH 1971 

C 

r ££J n k!r i: 2£,-££JL^ INIMUM ALONG THE SEAPCh VECTOR USING THE 
C C-CLDEN SECTION SEARCH METHQC 

C CHITON/ SHARE/ X(IOO), CEL(IOO), A ( 1GO, IOC ) , N , M , NN,NF1 

1 t b r 1 

CCNHCN /VALUE/ F , G , PO ,R S I GMA , RJ(200), RHC 
CCHKCN/CRST/ DcLX ( 100 ) , DELXGdOO), RHJIN, RATIO, EPSI, 
1THETA0, 

2 RSIG1, Gl, X 1 ( 100 ) , X2(100), >2(100), XR2(100), 

3XP 1 (100 ) ,PP1, 

A FggtPlf FI, RJK200), D0T7, PGRAC(lOO), CIAC-(IOO), 

5 PR EV3, AC ELX, NTCTR, NUNINI, NPHASE, NSA7IS 
AES(CUMMY)=DABS(DUMNY) 

K S V = 1 
N 40 5 = 1 
F 2 1 = P0 
I S 1 
CCT7=0» 

CC 10 J = 1 , N 

CCT7=D_TT + DELX( J ) *QEL XO ( J ) 

GC 70 40 
CC 20 1*1, N 
C EL X( I ) =-OELX( I ) 

CCN7INUE 
N^C4=0 
NMNN + l 

C MN IS NCW NUMB. OF P0IN7S AFTER MIN ACHIEVED 
N TCTR=N TCTR + 1 
CC 50 1 = 1, N 
50 X 2 ( I ) =X ( I l 

FX1=P0 
N4C 1 = 0 
N AC 1=N4G1 + 1 
CC 70 1*1, N 
X (I ) = X2 (I ) +DELX ( I ) 

CALL EVALU 

1 MEANS SATIS. OF CONSTRAINT NT.PREV. 2MEANS NOCHANGE 2 
MEANS VIOLATION 

IF FC1NT IS NOT FEASIBLE GIVE IT AN ARBITRARILY HIGH VALU 
GC TO (540,90,80), NSATIS 
F>2=10.E25 
FC=10.E25 
GC 70 100 
CONTINUE 
FX2=P0 

IF ( PX1-PX2 ) 100,100,150 
IF (N401-2) 130,110,110 
CC 120 1*1, N 



10 

20 

20 

40 



60 

70 

C 

C 

c 

£0 



90 



100 

110 

120 



SG FAR COMPUTEC 



X 1 ( I ) =X ( I ) 

F 1= PX2 
C-C TO 420 
C ONLY ONE PCINT 
120 CC 140 1*1, N 

140 X2(I)=X2(I) 

FREV2=P>1 
GC 70 1 €0 

150 CC 160 1 = 1, N 

X 2 ( I ) =X 2( I ) 

160 CELX(i) =1.61803399*CELX(I ) 
F F EV3 = P XI 
F ) 1=PX2 
GC TO 60 

C GCLDEN SECTION SEARCH METHOD. 

C E VECTOR C-CES TO XI ( I ) 

170 FC= 1.E2 6 , 

N 40 4 = N404 + 1 
180 DC 190 1*1, N 



190 xni)sx(i) 

F 1= F0 

CC 2C0 1=1, N 

X (I )=.38196601 S MX1( I)-X3( I) )+X3(I) 

200 X2(I)=X(I) 

CALL EVALU 

C-C TC (540,270,210), NSATIS 

210 IF (N4C4.LT. 30) GC TC 170 

211 CONTINUE 

C THERE IS NC REFERENCE TC 211, T HE AECVE STATEMENT IS A 
C DUMMY STATEMENT 

C— IT IS POSSIBLE NO FEASIBLE POINT EXIST, IF NCT TRY MOVING 
C ON CELXO 

C-- IF IT IS NOT POSSIBLE TC MOVE ON CELXO THEN WE MUST EE 
C AT A SOLUTION CF NLP PROBLEM. 

IF (N404.GT.100) GO TO 240 
220 CC 230 1 = 1, N 

IF <AES<ABS<X3<I )/Xl(I) )-l 9 )„GT. l.E-7) GC TC 170 
230 CONTINUE 
240 GC TC (250,260), N405 

250 N 4C 5= 2 

C — TRY TC MCV E CN GRADIENT. 

NTCTR=NTCTR-1 
MN=MN-1 
GC TC 20 

260 WRITE (6,580) 

CALI TIMEC 
CALL OUTPUT (1) 

CALL REJECT 
STOP 22042 
C 

270 CONTINUE 
N4C4=0 
FX1=P0 

CC 230 1 = 1 , N 

280 X (I ) =0.28196 60 1*( X 1 ( I )-X2(I))+X2(I) 

CALL EVALU 

GC TO (540,290,220), NSATIS 



290 P X 2 = P0 

N 40 1= 1 

300 N4C1=N401+1 

IF (N401-25) 340,210,310 
210 KSM2 

IF (N4C1-40) 320,460,460 
-pH r r " ^ 0 1 = 1, N 

IF (ABS < X 2 ( I ) / X ( I ) - 1 .0 ) .GE. l.E-7 ) GO TO 240 
220 CONTINUE 
GC TO 460 

340 IF (ABS (PX1/PX2-1. ).LE. l.E-7) GC TC 460 

C FRCM^LEFTGRIGFT* X3? it ( PRE V2 ) X2 ( 1 ) ( PX 1 )X ( 1 ) P X2 X1C)P1 
250 CC 260 1 = 1, N 
260 X 1 ( I ) =X ( I I 
C THROW AWAY RIGHT PART 



F 1 = PX2 

STS^^PlisHLSfxx.ii-xauM^m 

C TEMPORARILY IN X STORAGE 
CALL EVALU 

C-C TO ( 540,280, 170 ), NSATIS 
380 CONTINUE 

C SWITChSeCTCRS TO PROPER POSITION 
F X 1 = P0 

CC 290 1=1, N 
X X = X2 ( I ) 

X2 (I ) =X (I ) 

290 X(I)=XX 

GC TC 200 

C L p FT SICE TOSSED AWAY 
C — "CHANGES FOR NONUNI MCCAL FN 
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C-- GC 7 C ThRCW AWAY 



400 

410 

420 



420 

44C 



450 



C The 
460 



470 



480 

490 

500 

C IF 

510 

520 

520 

C ARE 
540 

550 



C 

c -• 



RIGI-T IN THIS CASEINIT VAL LT FIB FT 
350,250,410 



COMPUTE MI 



IF (PREV2-PX2) 

CC 420 1-1 »N 
X2(I)=X2(I> 

X2 ( I ) = X (I ) 

FFEV3=PX1 
F)1=PX2 
CC 440 1 = 1, N 

> ( I ) = 0. 28 156601 *( XI ( I l-X 2 ( I ) ) +X 2 ( I ) 

CALL EVALU 

C-C TC (540,450, 170), NSATIS 
CONTINUE 
FX2=PQ 
GC TC 300 

INTERIOR PGINTS NOW GIVE EQUAL VALUE FOR F, 

CC 470 1=1, N 
CELXOII ) = X ( I ) 

X (I ) = (OELXC( I ) + X2d) )*0.5 

CONTINUE 

CALL EVALU 

GC TO (460,490). KSW 

IF (ASS (PQ/PX1-1. ) »GT«l«E-7 ) GO TC 520 
GC TO (500,510), ISW 
IF (P0.LT.P21) GO TO 510 
I S W = 2 

P-FLNCTICN Cl ON, T GC DOWN TRY NEC- VECT. 

GC TO 20 
RETURN 

CC 530 1=1, N 
X ( I ) = DEL XO ( I ) 

C-C TO 2 50 

WE NOW IN FEASIBILITY PHASE 
CC 550 I = 1 » M 
IF ( R J ( I ) ) 560, 5oO, 550 
CCNTINLE 
NSATIS=4 
PETLRN 

PROBLEM HAS BECOME FEASIBLE 
F - FUNCTION CHANGES IF A CONSTRAINT BECOMES FEASIBLE 



560 


MN-C 






CC 570 


1 = 1, M 


570 


R„1 (I ) = 


PJ(I) 




RETURN 




£ 

580 


FCFM AT 


( 8QH 



OPT CAN-T FINC A 



1G IVES A 
ENC 



FEASIBLE POINT, THAT 
INCTICN. ) 
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GLOSSARY FOR CUBIC INTERPOLATION OPT SUBROUTINE 



ADELX 

AL 

BL 

COTES 

D 

DELX 

DELXO 

DLXl 

DLX2 

DOTS 

DOTT 

DOTTA 

DOTTB 

EPSO 

I 

J 

K 

KTER 

KTR 

M 

MN 



Magnitude of the search direction vector 

Lower step length 0 T 

L 

Upper step length 0^ 

Cosine of the angle between the gradient and search 
direction vectors 

Factor by which increment is multiplied 
Search direction vector s 

Gradient vector of the penalty function VP(X) 

Gradient vector of penalty function for lower step 
length VP(X3) 

Gradient vector of penalty function for upper step 
length VP(X2) 

Directional derivative of interpolated step length 
times ADELX P' (X) 

Directional derivative of initial starting point 
times ADELX P' (XI) 

Directional derivative of lower step length times 
ADELX P' (X3) 

Directional derivative of upper step length times 
ADELX P ' (X2) 

Stopping criteria tolerance for cosine test 
Index used in DO loops 
Index used in DO loops 
Index used in DO loops 

Number of points evaluated in bracketing the minimum. 

Number of cubic interpolations performed 

Number of inequality constraints 

Number of moves in search for solution of a 
subproblem 
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N 


Number of variables 


NRED 


Number of consecutive step length reductions 


NSATIS 


Indicates whether constraints are satisfied 


NTCTR 


Number of the point on which the program is working 


N405 


Switch showing that a feasible point on the 
direction vector could not be found and the 
negative was tried. 


PO 


Penalty function value at current xvector P(X) 


PX1 


Penalty function value at lower step length P(X3) 


PX2 


Penalty function value at upper step length P(X2) 


Q 


Coefficient used in computing interpolated step length 


RJ 


Vector of current values of the constraints 


RJ1 


Vector of previous values of the constraints 


SMAG 


Magnitude of the gradient of the penalty function 


TAL 


Current step length 0 


TINC 


Increment by which step's length increased or 
decreased 


X 


Current X vector 


XI 


Initial starting point X vector 


X2 


X vector of upper step length x q +0 u s 


X3 


X vector of lower step length X Q +0 L S 


Z 


Coefficient used in computing interpolated step 
length 
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A 

c 



„ _ S LE RCUT INE OPT 

r I^fui C £ T A?y^ C y T J NE USES CUBIC INTER PC LATION TC FIND THE 
C NIMPLN ALONG THE GIVEN SEARCH VECTOR 
IMPLICIT R 5AL*8 ( 4 -H Q*"Z) 

CCNMON/SHARE/ X(IOO) 1 , CEL(IOO), A ( IOO , 1 00 ) , N , M , NN ,N F 1 

1 * n y l 

£?EM,/)(H U?/ ^ F » G » PO t RSIGMA » R J ( 2G0 ) * R HO 
CCNNCN/CRsT/ DELX(IOO), DcLXO(lOO), RHGIN, RATIO, EFSI, 
1 i HETAO, 

2 FSIG1, Gl, XI ( 100 ) , X 2 ( ICO ) , X2<100), XR2(100), 
3XR1 (100 ) , PR1, 

i R J 1 ( 20C ) , DOTT, PGRAC(IOC), CIAG(IOO), 

5 n PREV3, ADfcLX, NTCTR , NUMINI, NPHASE, NSATIS 
CINENSICN CLX1( 100 ) ,CLX2( 100 ) 

A ES ( CUMMY ) =CABS ( DUNMY ) 

S CRT ( DUMMY J-OSQRT ( DUMMY ) 

E FSC=1» E-7 
CCTS=0.0 
NAC 5= 1 
GC TO 5 
- CC A 1=1, N 

CELX( I) =-DELX( I ) 

CONTINUE 
INITIALIZATION 
NMNN + 1 
NTCTR=NTCTR+1 
CCTTA=0.0 
TAL=0.0 
AL=C.O 
C = 2.0 
K T E P = 0 
NFEC=0 
TINOl.O 
S N AG=Oo C 
CC 10 K = 1 , N 
SNAG=SNAG+CELX(K )**2 
CCTTA=CCTTA-DELX(K)*DELXO(K) 

X 1 ( K ) =X ( K ) 

CLXl(K) =DELX0(K) 

X 2 ( K ) =X ( K ) 

CCTT=-CCTTA 
PX1=P0 

C INITIAL BRACKETING BEGINS 
25 T AL=TAL*T INC 
2C CC 25 K = 1 , N 
25 X(K)=X1(K)+DELX(K)*TAL 
KT ER = KT ER + 1 
CALL EVALU 

C-C TG ( 170,50 ,40) , NSATIS 
C RECLCE STEP SIZE 
AC C =0 • 5 

NREC=NR ED+1 
TINC=(T AL-AL )*D 
TAL=TAl-TINC 
C-C TC 50 

C SECCNC FEASIBLE POINT FCUND 
5C B l =TAL 
F X 2 = P0 

CALL GR AD ( 2 ) 

CCTTB=0.0 
CC 60 K = 1 , N 
X 2 ( K ) =X ( K ) 

C L X 2 ( K ) =DELX0 ( K ) „ , , 

6C CC7T6=GCTTfi-DEL X ( K ) *D EL XO < K ) 

C CHECK STOPPING CRITERIA OR HAYBE REVERSE SEARCH DIRECTION 
I F (ADELX.EC.O .0) GC TO 91 
I F (SMAG .EC. 0.0 ) GC TO 91 
CCT ES= ABS ( DOTTB )/(SQRT<SMAG)*ACELX) 

IF (COT ES.LTaEPSQ ) GO TO 91 
IF (KTER.GT.100) GO TC 280 
I F (NREC .GT.20) GO TO <c80 
CC 75 K = 1 , N 



1C 



67 



75 CCIvfiNU£ BS<X1<K>/)<2<K * >_1 * * - GT - 1-E-7J GO TG 76 
GC TO 260 

76IF(£8$(A6S(PX2/PX1 )-l. ) .L E. l.E-7 ) GQ TO 91 
C IF CEPIVATIVE IS POSITIVE THEN THE MINIMUM IS BFACKETEC 
C AND CLEIC INTERPOLATION CAN BE PERFORMED, IF NC7 MAKE THIS 
C PCINT TFE LOWER STEP LENGTH AND STEF CUT FARTHER 
IF(COTTE.GT.O.O) GO TO 100 
N P E C = 0 
T INC = T I NC*D 
CC 60 K*1,N 
DLX 1( K) =DLX2(K ) 

EC X 2 ( K ) =X 2 ( K ) 

F > 1 = PX2 
AL = EL 

CC7 TA = CCTTB 
GC TO 25 

C OF THE TWO FEASIBLE POINTS RETURN WITH THE SMALLEST 

51 IF(FX2.LE.PX1) GC TO 220 
PC=PX1 

CC S 2 K = 1 , N 
X (K ) = X3 (K ) 

CELXOtK )=DLX1(K) 

5 2 CONTINUE 
GC TO 220 

C CLEIC INTERPOLATION FOLLOWS **************** *** 4*4** * * **** 
ICC K7F=0 
1C 5 CONTINUE 

Z=2.*( ( FX1-PX2 ) / ( EL-AL) ) +DOTT A + CCTTB 
C=SCRT(Z**2-DOTTA*CCTT8 ) 

TAL = 8L-(BL-AL)=M DCTTB + Q-Z ) / ( D07T E-CCTTA +2 • ) 

CC 130 K»1,N 

12C X <K) = X1 <K) +DELX(K)=MTAL 
CALL EVALU 

GC TO ( 170,135,40) , NSATIS 
125 CALL GR AD ( 2 ) 

CC7S=0.G 
CC 140 K = 1 , N 

CCTS=CCTS-CELX(K)*OELXO(K) 

1 AC CONTINUE 

C CHECK STOPPING CRITERION GR IF SEARCH DIRECTION SHOULD BE 
C REVERSEC 

I F(ADELX.EC.O.O) GO TO 471 
IF(SMAC-.EC.O.O) GC TO 471 
CCTE$=AES( CCTS)/ (SCRT(SMAG) *ADELX ) 

IF (CGTES.LT. EPSO ) GO TG 471 

IF(ABS(ABS (PX2/P0 ) -1. ) .L E. l.E-7 ) C-0 TO 471 

KTF=KTR+1 

IF (KTR.GT.10) GO TO 471 
CC 143 K=1,N 

IF(AE$(A8S(X2(K)/X(K) ) -1 . ) .GT. 1. E-7 ) GO TC 163 

142 CONTINUE 
GC TC 220 

162 IF (COTS .C-T. 0.0) GO TO 150 
C THROW AWAY LEFT SIDE 
A L = TAL 
PX1=PC 
CCTTA=CCTS 
CC 145 K = 1 ,N 
CLX1(K)=DELX0(K) 

145 X 2 ( K ) =X (K ) 

GC TO 1 05 

C THFCW AWAY RIGHT SIDE 
1 5 C 6 L = TAL 
FX2=P0 
CCTTe=CCTS 
CC 160 K= 1 »N , 

CLX2(K)=DELX0(K) 

16C X 2 ( K ) *X ( K ) 

GC TO 105 

170 CC 180 K = 1 ,M ^ 

I F (RJ (K ) ) 190,190,160 
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1 EG CONTINUE 

N SAT I S = 4 
RETURN 
ISC NN = C 

CC 200 K = 1,M 
2C0 RJ1(K)=RJ(K) 

RETURN 

RETURN WITH SMALLEST VALUE 

471 I F ( F0.LE.PX1 ) GO TO 474 
IF(FX1.LE.PX2) GO TO 475 

472 FC = FX2 

CC 473 K = i,N 
X (K ) = X2 (K ) 

C ELXO ( K ) = DUX2( K ) 

472 CCNTINUE 
GC TO 220 

474 IF1F0.LE.PX2) GO TC 220 
GC TC 472 

475 PC=PX1 

CC 476 K = 1,N 
X (K )=X3 (K) 

CELXO(K )=DUX1(K) 

476 CCNTINUE 
22C CCNTINUE 
225 CCNTINUE 

RETURN 

REVERSE SEARCH DIRECTIONS 

2 EC GC TO ( 290,300) ,N405 
29C N 40 5=2 

NTCTR=NTCTR-1 
NN=NN-l 
GC TO 3 

3CC fcFITE(6,580) 

CALL TI NEC 
CALL CUTPU'T(l) 

CALL REJECT 

5 £0 FORMAT ( 80H OPT CAN-T FIND A FEASIBLE PC IN T , THAT 
1C-IVES A L G to E R VALUE UF THE F-FLNCTIQN. ) 

STOP 22G42 
ENC 
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GLOSSARY FOR QUADRATIC INTERPOLATION OPT SUBROUTINE 



AL 


Interpolated step length 


D 


Factor by which increment is multiplied 


DELX 


Vector indicating the direction of move in one 
dimensional optimization S 


DELXO 


Gradient vector of the penalty function VP 


DOTT 


Directional derivative of initial starting point 
times search direction magnitude P'(X1) 


I 


Index used in DO loops 


J 


Index used in DO loops 


K 


Index used in DO loops 


KTER 


Number of points evaluated in bracketing the 
minimum 


KTR 


Number of quadratic interpolations performed 


M 


Number of inequality constraints 


MN 


Number of moves in search for solution of a 
subproblem 


N 


Number of variables 


NFIRS 


Indicates if more than two feasible step lengths 
found 


NRED 


Number of consecutive step length reductions 


NSATIS 


Indicates whether constraints are satisfied 


NTCTR 


Number of the point on which the program is working 


N405 


Switch showing that a feasible point on the direction 
vector could not be found and the negative was tried 


PO 


Current penalty function value P(X) 


PX1 


Penalty function value of lower step length P(X3) 
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PX2 

RJ 

TAL 

TALI 

TAL 2 

TINC 

X 

XX 

XI 
X2 
X3 



Penalty function value of upper step length 

Vector of previous values of the constraints 

Current step length 0 

Lower step length 0 T 

Li 

Upper step length 0 u 

Increment by which step length increased or 
decreased 

Current X vector Xq+ 0s 

Temporary storage used for switching values 
Initial starting point vector Xq 
X vector at upper step length x q + ®u S 
X vector of lower step length x q + ®L S 
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SIEROUTINE OPT 

C TUS CPT SUBROUTINE FINDS THE MINIMUM ALONG A GIVEN SEARCH 
C VECTOR USING CUADRATIC INTERPOLATION 
IMPLICIT SEALES ( A-H «0-Z ) 

COMMON/ SHARE/ X(IOO), DEL(IOO), A { 100 , 100 ) , N , M , MN,NF1 
1,NM1 

COMMON /VALUE/ F , G , PO , RS I C-M A , RJ(20Q), RFC 
COMMON/ CRST / DELX(IOO), DELXG(IOO), RHQIN, PATIO, EPSI, 
1TFETA0, 

2 PSIC-1, Gl, XI ( 100 ) , X2 ( 100 ) , >3(10'')), X02(iOC), 
2XPK100 ),PRl, 

A F F 2 , P 1 , FI, RJK200), DOTT , PGRAC(IOO), DIAG(IOO), 

5 PPEV3 , AOELX, NTCTR , NUMINI, NPHASE , NSATIS 
AES (CUMMY)=DABS( DUMMY) 

N AO 5= 1 
CCTT=C. 

C C 2 K “ 1 N 

2 CCTT=CCTT+CELX(K)*CELXO{K) 

C-C TO 5 
3 DC A 1 = 1, N 

C ELX ( I ) =-DELX( I ) 

CONTINUE 
C INITIALIZATION 
M N = MN + 1 
NTCTR = N TCTR + l 
T AL=0*0 
T AL 1 = 0*0 
C = 2 • 

KTER=0 
NREC=0 
T INC= 1. 0 
DC 7 K= 1, N 
X 1 ( K ) =X ( K ) 

X 2 ( K ) =X (K) 

7 F X 1 = P0 

N F I RS =0 

C START BRACKETING PROCEDURE 
T AL = TAL+T INC 
10 DC 13 K = 1 » N 
12 X (K ) = X1 (K)+DELX(K )*TAL 
KTSR=KTcR+l 
CALL EVALU 

GO TO ( 170,20,20 ) , NSATIS 
C RECUCE STEP SIZE 

20 C = C • 5 

IF (NREC.GT.20) GO TO 280 

NR EC = NR ED + 1 

TINC = (T AL-TAL1)*D 

T AL = TAL-T INC 

C-C TO 1C 

C HAVE 2 CR MORE FEASIBLE POINTS BEEN FOUND 
30 I F ( NF IP S • EC* 1 ) GO TC AO 
N F I RS = 1 

21 CC 22 K = 1 * N 

22 X 2 ( K ) =X ( K i 
F X 2 = PC 

IF (PX2.C-E.PX1) GO TO 20 
C INCREASE STEP SIZE 
35 T INC=TI NC*D 

T AL =TAL+T INC 
N R E C = 0 
C-C TO 10 
AO CONTINUE 

C RECRCER PCINTS „ _ 

IF(TAL.LT.TAL2) GO TO 5 A 

I F ( FO.L T.PX2 ) GO TO 50 

XX=TAL 

7 AL=TAL2 

T A L 2 = XX 

XX = FO 

PC = PX2 
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F >2=XX 

CC 45 K = 1 , N 

XX=X(K) 

X ( K ) = X2 ( K ) 


45 


X 2 ( K ) =X X 
GC TO 54 


46 


TAL1=TAL 
PX1=P0 
CC 47 K = 1 , N 
X 3 ( K ) = X ( K ) 

T AL=TAL 2 
GC TO 25 


47 


50 


T AL 1 = T A L2 
F > 1 = PX2 
CC 51 K*1,N 
X 3 ( K ) =X 2 ( K ) 


51 


X 2 ( K ) =X ( K ) 
T AL2=TAL 
FX2=P0 



52 

c -a 

c 

c 

54 

c c 

56 



C CI-ECK STOPPING CRITERION OR IF SEARCH OIRECTION SHOULC eE 
C REVERSEC 

IF(AeS(AB$(PO/PXl )-!.). LT.l.E-7) GC TO 220 
CC 52 K = 1 » N 

IF( AB$( ABS(X1(K)/X(K) )-l. ).GT.l.E-7) GO TC 55 
CCNTINLE 
GC TO 2 60 
CONTINUE 
GC TO 35 

CHECK STOPPING CRITERION OR IF SEARCH DIRECTION SHOULC eS 
REVERSEC 

IF(AeS(ABS(PO/PXl)-l.). LT.l.E-7) C-0 TO 220 
CC 55 K = 1 , N 

IF (ABS (ABS(XKK) / X ( K i )-I • ).GT.I.E-7) GO TC 56 
CCNTINUE 
GC TO 2E0 
CCNTINLE 

IF (KTER.GE.IOO) GO TO 280 
IF(P0.GT.PX2) GO TO 46 
IF(PO.GE.PXl) GO TO 31 
CC 57 K = 1 , N 

IF(ABS(A8S(X2(K)/X(K> ) -1. ) . G 6. 1 . E-7 ) GO TC 60 
57 CCNTINUE 
GC TO 220 

C CLACRA7IC INTERPOLATION BEGINS*************************** 

60 KT R =0 

61 k!P=KTS+1 

AL=C. 5* ( ( TAL**2-TAL2**2)*PX1+ (TA L 2*^2— T AL 1* *2 ) *P0 + 

1 (TAL1**2-TAL**2 

2 ) *FX2)/ ( (TAL-TAL 2 )*PXl+ (T AL 2-TA L 1 ) *P0 + ( T AL 1 -T AU * P X 2 ) 
IF(TAL.C-T.AL) GO TC 75 

C THROW AWAY LEFT SIDE 
T AL 1 = TA L 
FX1=P0 
CC 71 K = 1 » N 
71 X3(K)=X(K) 

GC TO 78 

C TI-PCW A ^AY FIGHT sice 
75 TAL2=TAL 

P X 2 = P0 
CC 77 K = 1 , N 

77 X 2 ( K ) =X ( K ) 

78 TAL=AL 

CC 79 K = 1»N 

79 X (K ) = X 1 (K)+DELX(K )*TAL 
CALL EVALU 

GC TO (170,80,85), NSATTS 

C C CHECK C STOPPINd N CRITERICN OR IF SEARCH DIRECTION SHOULD 85 
c rfC=rcec 

IF(AES(ABS(X2(K)/X(K) )-l. J.GE.l.E-7) GO TC 83 
82 CCNTINUE 

GC TG 220 
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82 It !5SS ( £S S <P0/PX1 )-l„).LT.l.c-7) GOTO 220 

IMKTR.GT.20) GQ TG 220 
„„ GC TO 61 

6 5 N R E C = 0 

KT £R=0 
GC TO 20 

17C CC 180 K = 1,N 

IF(RJ(K)» 190,190,180 
160 CONTINUE 
N S AT I S= 4 
RETURN 
190 NN=0 

CC 200 K = 1,M 
2CC RJ1(K)=PJ(K) 

RETURN 

220 CONTINUE 

C RETURN WITH SMALLEST VALUE 
IF(PXl.GE.PO) GQ TC 240 
PC=PX1 

CC 230 K=1,N 
230 X (K )=X3 (K ) 

240 CONTINUE 
P ET UR N 

C REVERSE SEARCH DIRECTIONS 
2 EC GC TO ( 290,300) ,N405 
290 N40 5 = 2 

NTCTR=NTCTR-1 
NN=NN-1 
GC TG 3 

300 W P I TE ( 6 , 580 ) 

CALL TI NEC 
CALL OUTPUT ( 1 ) 

CALL REJECT 

5 EC FC FNAT ( 80H OPT CAN-T FIND A FEASIBLE POINT, ThAT 
1GIVES A LOWER VALUE OF THE P-FUNCTIGN. ) 

STOP 22C42 
ENC 
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